#### setting environment ####
# install.packages("estimatr")
require(estimatr)

# function to draw the point estimate and the confidence interval of a group mean
point.draw <- function(x, y, pch, bg = NA) {
  res <- lm_robust(y ~ 1)
  segments(x, res$conf.low, x, res$conf.high)
  points(x, res$coefficients, pch = pch, bg = bg)
}

# function to test the significance of the sum of two coefficients
sum.coef.test <- function(result, var1, var2) {
  est <- result$coefficients[var1] + result$coefficients[var2]
  se <- sqrt(result$vcov[var1, var1] + result$vcov[var2, var2] + 
               2 * result$vcov[var1, var2])
  t.stat <- est / se
  p.val <- 2 * (1 - pt(abs(t.stat), result$df.residual))
  out <- c(est, se, t.stat, p.val)
  names(out) <- c("Estimate", "Std. Error", "t value", "Pr(>|t|)")
  out
}

# function to test the significance of the difference between two coefficients
dif.coef.test <- function(result, var1, var2) {
  est <- result$coefficients[var1] - result$coefficients[var2]
  se <- sqrt(result$vcov[var1, var1] + result$vcov[var2, var2] - 
               2 * result$vcov[var1, var2])
  t.stat <- est / se
  p.val <- 2 * (1 - pt(abs(t.stat), result$df.residual))
  out <- c(est, se, t.stat, p.val)
  names(out) <- c("Estimate", "Std. Error", "t value", "Pr(>|t|)")
  out
}

# function to compute a pooled standard deviation
pooled.sd <- function(x1, x2) {
  n1 <- length(x1)
  n2 <- length(x2)
  sqrt(((n1 - 1) * var(x1, na.rm = TRUE) + 
          (n2 - 1) * var(x2, na.rm = TRUE)) / (n1 + n2 - 2))
}

#### data ####
raw.data <- read.csv("replication_data.csv")

# number of respondents in the original data
nrow(raw.data)

# exclude inattentive respondents and
# respondents who lacked information required for blocked
data <- subset(raw.data, block != "Not target of blocked randomization")

# number of respondents in the cleaned data
nrow(data)

#### factual manipulation checks ####
# Online Appendix D, Table A.1(a)
table(data$treatment, data$MC1)
round(prop.table(table(data$treatment, data$MC1), margin = 1), 2)
table(data$treatment)
table(data$MC1)
round(prop.table(table(data$MC1)), 2)
sum(table(data$MC1))

# Online Appendix D, Table A.1(b)
table(data$treatment, data$MC2)
round(prop.table(table(data$treatment, data$MC2), margin = 1), 2)
table(data$treatment)
table(data$MC2)
round(prop.table(table(data$MC2)), 2)
sum(table(data$MC2))

#### average outcome variables in each experimental group ####
## Figure 2
cairo_pdf("Figure_2.pdf", width = 6, height = 4, pointsize = 9, family = "Helvetica")
layout(matrix(1:6, 2, 3, byrow = TRUE))
par(mar = c(3, 2.5, 4, 0.5), lwd = 0.5)
plot(NULL, NULL, type = "n", bty = "n", xlim = c(0.5, 4.5), ylim = c(3, 4), 
     main = "Consensus mobilization\n", xlab = "", ylab = "", xaxt = "n", yaxt = "n")
abline(h = seq(3, 4, 0.25), lty = 3, col = "gray")
point.draw(1, data$consensus[data$treatment == "Control"], 4)
point.draw(1.85, data$consensus[data$treatment == "Treatment 1"], 19)
point.draw(2.15, data$consensus[data$treatment == "Treatment 2"], 21, bg = "white")
point.draw(2.85, data$consensus[data$treatment == "Treatment 3"], 19)
point.draw(3.15, data$consensus[data$treatment == "Treatment 4"], 21, bg = "white")
point.draw(3.85, data$consensus[data$treatment == "Treatment 5"], 19)
point.draw(4.15, data$consensus[data$treatment == "Treatment 6"], 21, bg = "white")
mtext("Control\n", side = 1, at = 1, line = 1, cex = 0.7)
mtext("Personal\n", side = 1, at = 2, line = 1, cex = 0.7)
mtext("Informa-\ntional", side = 1, at = 3, line = 1, cex = 0.7)
mtext("Motiva-\ntional", side = 1, at = 4, line = 1, cex = 0.7)
axis(2, at = seq(3, 4, 0.25), lwd = 0.5, cex = 0.7)
plot(NULL, NULL, type = "n", bty = "n", xlim = c(0.5, 4.5), ylim = c(25, 55), 
     main = "Action mobilization\n(Donation)", 
     xlab = "", ylab = "", xaxt = "n", yaxt = "n")
abline(h = seq(25, 55, 10), lty = 3, col = "gray")
point.draw(1, data$donation[data$treatment == "Control"], 4)
point.draw(1.85, data$donation[data$treatment == "Treatment 1"], 19)
point.draw(2.15, data$donation[data$treatment == "Treatment 2"], 21, bg = "white")
point.draw(2.85, data$donation[data$treatment == "Treatment 3"], 19)
point.draw(3.15, data$donation[data$treatment == "Treatment 4"], 21, bg = "white")
point.draw(3.85, data$donation[data$treatment == "Treatment 5"], 19)
point.draw(4.15, data$donation[data$treatment == "Treatment 6"], 21, bg = "white")
mtext("Control\n", side = 1, at = 1, line = 1, cex = 0.7)
mtext("Personal\n", side = 1, at = 2, line = 1, cex = 0.7)
mtext("Informa-\ntional", side = 1, at = 3, line = 1, cex = 0.7)
mtext("Motiva-\ntional", side = 1, at = 4, line = 1, cex = 0.7)
axis(2, at = seq(25, 55, 10), lwd = 0.5)
plot(NULL, NULL, type = "n", bty = "n", xlim = c(0.5, 4.5), ylim = c(2.75, 3.75), 
     main = "Action mobilization\n(Social media sharing)", 
     xlab = "", ylab = "", xaxt = "n", yaxt = "n")
abline(h = seq(2.75, 3.75, 0.25), lty = 3, col = "gray")
text(1, 3.25, "Not\napplicable")
point.draw(1.85, data$sharing[data$treatment == "Treatment 1"], 19)
point.draw(2.15, data$sharing[data$treatment == "Treatment 2"], 21, bg = "white")
point.draw(2.85, data$sharing[data$treatment == "Treatment 3"], 19)
point.draw(3.15, data$sharing[data$treatment == "Treatment 4"], 21, bg = "white")
point.draw(3.85, data$sharing[data$treatment == "Treatment 5"], 19)
point.draw(4.15, data$sharing[data$treatment == "Treatment 6"], 21, bg = "white")
mtext("Control\n", side = 1, at = 1, line = 1, cex = 0.7)
mtext("Personal\n", side = 1, at = 2, line = 1, cex = 0.7)
mtext("Informa-\ntional", side = 1, at = 3, line = 1, cex = 0.7)
mtext("Motiva-\ntional", side = 1, at = 4, line = 1, cex = 0.7)
axis(2, at = seq(2.75, 3.75, 0.25), lwd = 0.5)
plot(NULL, NULL, type = "n", bty = "n", xlim = c(0.5, 4.5), ylim = c(3, 4.5), 
     main = "Empathy\n", xlab = "", ylab = "", xaxt = "n", yaxt = "n")
abline(h = seq(3, 4.5, 0.5), lty = 3, col = "gray")
point.draw(1, data$empathy[data$treatment == "Control"], 4)
point.draw(1.85, data$empathy[data$treatment == "Treatment 1"], 19)
point.draw(2.15, data$empathy[data$treatment == "Treatment 2"], 21, bg = "white")
point.draw(2.85, data$empathy[data$treatment == "Treatment 3"], 19)
point.draw(3.15, data$empathy[data$treatment == "Treatment 4"], 21, bg = "white")
point.draw(3.85, data$empathy[data$treatment == "Treatment 5"], 19)
point.draw(4.15, data$empathy[data$treatment == "Treatment 6"], 21, bg = "white")
mtext("Control\n", side = 1, at = 1, line = 1, cex = 0.7)
mtext("Personal\n", side = 1, at = 2, line = 1, cex = 0.7)
mtext("Informa-\ntional", side = 1, at = 3, line = 1, cex = 0.7)
mtext("Motiva-\ntional", side = 1, at = 4, line = 1, cex = 0.7)
axis(2, at = seq(3, 4.5, 0.5), lwd = 0.5)
plot(NULL, NULL, type = "n", bty = "n", xlim = c(0.5, 4.5), ylim = c(3, 4.5), 
     main = "Anger\n", xlab = "", ylab = "", xaxt = "n", yaxt = "n")
abline(h = seq(3, 4.5, 0.5), lty = 3, col = "gray")
point.draw(1, data$anger[data$treatment == "Control"], 4)
point.draw(1.85, data$anger[data$treatment == "Treatment 1"], 19)
point.draw(2.15, data$anger[data$treatment == "Treatment 2"], 21, bg = "white")
point.draw(2.85, data$anger[data$treatment == "Treatment 3"], 19)
point.draw(3.15, data$anger[data$treatment == "Treatment 4"], 21, bg = "white")
point.draw(3.85, data$anger[data$treatment == "Treatment 5"], 19)
point.draw(4.15, data$anger[data$treatment == "Treatment 6"], 21, bg = "white")
mtext("Control\n", side = 1, at = 1, line = 1, cex = 0.7)
mtext("Personal\n", side = 1, at = 2, line = 1, cex = 0.7)
mtext("Informa-\ntional", side = 1, at = 3, line = 1, cex = 0.7)
mtext("Motiva-\ntional", side = 1, at = 4, line = 1, cex = 0.7)
axis(2, at = seq(3, 4.5, 0.5), lwd = 0.5)
plot(NULL, NULL, type = "n", bty = "n", xlim = c(0, 1), ylim = c(0, 1), 
     xlab = "", ylab = "", xaxt = "n", yaxt = "n")
legend("center", 
       legend = c("Control", "w/o justice frame", "w/ justice frame"), 
       pch = c(4, 19, 21), bg = c(NA, NA, "white"), bty = "n", cex = 1.2)
dev.off()

#### main results for consensus mobilization ####
# estimation of Model (1)
outcome1.test <- lm_robust(consensus ~ personal + informational + motivational + 
                             (personal + informational + motivational) : justice, 
                           data = data, fixed_effects = ~ block)
outcome1.test.summary <- summary(outcome1.test)

# number and percentage of missing outcomes
nrow(data) - outcome1.test.summary$nobs
round((nrow(data) - outcome1.test.summary$nobs) / nrow(data), 3)

# Table 2, columns 1-3
round(outcome1.test.summary$coefficients, 3)
round(dif.coef.test(outcome1.test, 1, 2), 3)  # beta_1 - beta_2
round(dif.coef.test(outcome1.test, 1, 3), 3)  # beta_1 - beta_3
round(dif.coef.test(outcome1.test, 4, 5), 3)  # delta_1 - delta_2
round(dif.coef.test(outcome1.test, 4, 6), 3)  # delta_1 - delta_3
outcome1.test.summary$nobs  # N. obs.

# effect size of the personal frame
round(outcome1.test$coefficients[1] / 
        pooled.sd(data$consensus[data$treatment == "Control"], 
                  data$consensus[data$treatment == "Treatment 1"]), 3)

# treatment effect of the motivational frame with the justice frame (beta_3 + delta_3)
round(sum.coef.test(outcome1.test, 3, 6), 3)

# treatment effect of the personal frame with the justice frame
round(sum.coef.test(outcome1.test, 1, 4), 3)

# effect size of the personal frame with the justice frame
round((outcome1.test$coefficients[1] + outcome1.test$coefficients[4]) / 
        pooled.sd(data$anger[data$treatment == "Control"], 
                  data$anger[data$treatment == "Treatment 4"]), 3)

#### main results for action mobilization: donation ####
# estimation of Model (1)
outcome2.test <- lm_robust(donation ~ personal + informational + motivational + 
                             (personal + informational + motivational) : justice, 
                           data = data, fixed_effects = ~ block)
outcome2.test.summary <- summary(outcome2.test)

# Table 2, columns 4-6
round(outcome2.test.summary$coefficients, 3)
round(dif.coef.test(outcome2.test, 1, 2), 3)  # beta_1 - beta_2
round(dif.coef.test(outcome2.test, 1, 3), 3)  # beta_1 - beta_3
round(dif.coef.test(outcome2.test, 4, 5), 3)  # delta_1 - delta_2
round(dif.coef.test(outcome2.test, 4, 6), 3)  # delta_1 - delta_3
outcome2.test.summary$nobs  # N. obs.

#### main results for action mobilization: social-media sharing ####
# estimation of Model (2)
outcome3.test <- lm_robust(sharing ~ (informational + motivational) * justice, 
                           data = data, fixed_effects = ~ block)
outcome3.test.summary <- summary(outcome3.test)

# number and percentage of missing outcomes
sum(data$treatment != "Control") - outcome3.test.summary$nobs
round((sum(data$treatment != "Control") - outcome3.test.summary$nobs) / 
        sum(data$treatment != "Control"), 3)

# Table 2, columns 7-9
round(outcome3.test.summary$coefficients, 3)
outcome3.test.summary$nobs

#### main results for emotional response: empathy ####
# estimation of Model (1)
outcome4.test <- lm_robust(empathy ~ personal + informational + motivational + 
                             (personal + informational + motivational) : justice, 
                           data = data, fixed_effects = ~ block)
outcome4.test.summary <- summary(outcome4.test)

# Table 3, columns 1-3
round(outcome4.test.summary$coefficients, 3)
round(dif.coef.test(outcome4.test, 1, 2), 3)  # beta_1 - beta_2
round(dif.coef.test(outcome4.test, 1, 3), 3)  # beta_1 - beta_3
round(dif.coef.test(outcome4.test, 4, 5), 3)  # delta_1 - delta_2
round(dif.coef.test(outcome4.test, 4, 6), 3)  # delta_1 - delta_3
outcome4.test.summary$nobs  # N. obs.

# treatment effect of the personal frame with the justice frame
round(sum.coef.test(outcome4.test, 1, 4), 3)

# effect size of the personal frame with the justice frame
round((outcome4.test$coefficients[1] + outcome4.test$coefficients[4]) / 
        pooled.sd(data$empathy[data$treatment == "Control"], 
                  data$empathy[data$treatment == "Treatment 4"]), 3)

#### main results for emotional response: anger ####
# estimation of Model (1)
outcome5.test <- lm_robust(anger ~ personal + informational + motivational + 
                             (personal + informational + motivational) : justice, 
                           data = data, fixed_effects = ~ block)
outcome5.test.summary <- summary(outcome5.test)

# Table 3, columns 4-6
round(outcome5.test.summary$coefficients, 3)
round(dif.coef.test(outcome5.test, 1, 2), 3)  # beta_1 - beta_2
round(dif.coef.test(outcome5.test, 1, 3), 3)  # beta_1 - beta_3
round(dif.coef.test(outcome5.test, 4, 5), 3)  # delta_1 - delta_2
round(dif.coef.test(outcome5.test, 4, 6), 3)  # delta_1 - delta_3
outcome5.test.summary$nobs  # N. obs.

# effect size of the personal frame
round(outcome5.test$coefficients[1] /
        pooled.sd(data$anger[data$treatment == "Control"],
                  data$anger[data$treatment == "Treatment 1"]), 3)

# treatment effect of the personal frame with the justice frame
round(sum.coef.test(outcome5.test, 1, 4), 3)

# effect size of the personal frame with the justice frame
round((outcome5.test$coefficients[1] + outcome5.test$coefficients[4]) / 
        pooled.sd(data$anger[data$treatment == "Control"], 
                  data$anger[data$treatment == "Treatment 4"]), 3)

#### additional tests for the treatment effect of the justice frame ####
## consensus mobilization
# estimation of Model (3)
outcome1.add.test <- lm_robust(consensus ~ informational + motivational + justice, 
                               data = data, subset = treatment != "Control", 
                               fixed_effects = ~ block)
outcome1.add.test.summary <- summary(outcome1.add.test)

# Online Appendix E, Table A.2, columns 1-3
round(outcome1.add.test.summary$coefficients, 3)
outcome1.add.test.summary$nobs  # N. obs.

## action mobilization: donation
# estimation of Model (3)
outcome2.add.test <- lm_robust(donation ~ informational + motivational + justice, 
                               data = data, subset = treatment != "Control", 
                               fixed_effects = ~ block)
outcome2.add.test.summary <- summary(outcome2.add.test)

# Online Appendix E, Table A.2, columns 4-6
round(outcome2.add.test.summary$coefficients, 3)
outcome2.add.test.summary$nobs  # N. obs.

## action mobilization: social-media sharing
# estimation of Model (3)
outcome3.add.test <- lm_robust(sharing ~ informational + motivational + justice, 
                               data = data, fixed_effects = ~ block)
outcome3.add.test.summary <- summary(outcome3.add.test)

# Online Appendix E, Table A.2, columns 7-9
round(outcome3.add.test.summary$coefficients, 3)
outcome3.add.test.summary$nobs  # N. obs.

# effect size of the justice frame
round(outcome3.add.test$coefficients[3] / 
        pooled.sd(data$consensus[data$treatment != "Control" & data$justice == 0], 
                  data$consensus[data$treatment != "Control" & data$justice == 1]), 3)

## emotional response: empathy
# estimation of Model (3)
outcome4.add.test <- lm_robust(empathy ~ informational + motivational + justice, 
                               data = data, subset = treatment != "Control", 
                               fixed_effects = ~ block)
outcome4.add.test.summary <- summary(outcome4.add.test)

# Online Appendix E, Table A.3, columns 1-3
round(outcome4.add.test.summary$coefficients, 3)
outcome4.add.test.summary$nobs  # N. obs.

## emotional response: anger
# estimation of Model (3)
outcome5.add.test <- lm_robust(anger ~ informational + motivational + justice, 
                               data = data, subset = treatment != "Control", 
                               fixed_effects = ~ block)
outcome5.add.test.summary <- summary(outcome5.add.test)

# Online Appendix E, Table A.3, columns 4-6
round(outcome5.add.test.summary$coefficients, 3)
outcome5.add.test.summary$nobs  # N. obs.